{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "from sklearn.gaussian_process import GaussianProcessRegressor\n",
    "from sklearn.gaussian_process.kernels import RBF, ConstantKernel as C\n",
    "from sklearn.gaussian_process.kernels import Matern"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.array([-2, 0, 4, 2]).reshape(-1,1)\n",
    "y = np.array([1, 0.3, -1, 0]).reshape(-1,1)\n",
    "\n",
    "x_test = np.linspace(-5, 5, 100).reshape(-1,1)\n",
    "\n",
    "kernel = RBF(1, (1e-2, 1e2))\n",
    "gp = GaussianProcessRegressor(kernel=kernel, n_restarts_optimizer=9)\n",
    "\n",
    "samples = gp.sample_y(x_test, 3)\n",
    "\n",
    "#gp.fit(x, y)\n",
    "#opt_params = gp.get_params()\n",
    "#print \"Opt param\", opt_params\n",
    "\n",
    "#log_lik = gp.log_marginal_likelihood()\n",
    "\n",
    "# Make the prediction on the meshed x-axis (ask for MSE as well)\n",
    "#y_pred, sigma = gp.predict(x_test, return_std=True)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD8CAYAAABjAo9vAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzs3XdcVfX/wPHXYW9kCIiIE3HgRnHPXJmaZubeWs5KzWyq9c3KsrLUzK25R6ZpztwLBZQl4kKGgICoyIZ7z++PU/2sUNa93Hvh83w8fFByzzlvFd587me835IsywiCIAjlh5GuAxAEQRA0SyR2QRCEckYkdkEQhHJGJHZBEIRyRiR2QRCEckYkdkEQhHJGJHZBEIRyRiR2QRCEckYkdkEQhHLGRBcPdXZ2lmvUqKGLRwuCIBiswMDAFFmWKxf2Op0k9ho1ahAQEKCLRwuCIBgsSZKii/I6MRUjCIJQzojELgiCUM6UOrFLklRNkqQTkiRFSJIULknSm5oITBAEQSgZTcyx5wOzZFkOkiTJFgiUJOmoLMvXNHBvQRAEoZhKPWKXZTlBluWgP//7CRABVC3tfQVBEISS0egcuyRJNYBmgL8m7ysIgiAUncYSuyRJNsBu4C1ZltMK+PwkSZICJEkKSE5O1tRjBUEQhH/RyD52SZJMUZL6ZlmWfynoNbIsrwRWAvj6+op+fIIglEvZeSoCox+Skp7Dg/RcsvJUeLva0tSzEs425mUSQ6kTuyRJErAGiJBl+ZvShyQIgmB4btx/whb/GH4JiiMtO7/A13g4WLJoUGPa1nbWaiyaGLG3A0YCoZIkXf3z996XZfl3DdxbEARBr2Xk5DNvXzi7AuMwMzaip48bA5tXxdPRCkcrM8xMjLiWkMaVmIdcjX2Eq52F1mMqdWKXZfksIGkgFkEQBIMSdu8x07de4e6DDCZ3rs2E9jVxKmC6pWUNR1rWcCyzuHRSK0YQBMHQ/RIUx7u7Q3CyNmfLhNa0qe2k65D+JhK7IAhCMR0ISWD2zmBa13Ji2bDmOFib6TqkfxCJXRAEoRhORCbx1vYrNPd0YPVoX6zM9C+NiiJggiAIRXT5bipv/ByIt5sta8e21MukDiKxC4IgFMmD9BymbA6iaiVLNoxthZ2Fqa5Deib9/HEjCIKgR2RZZvbOYB5n5bFhbKsCd77oEzFiFwRBKMS6c3c5EZnM+73r0cDdTtfhFEokdkEQhOcIj3/MFwev062eC6Pb1tB1OEUiErsgCMIzqNQy7+wMwd7KlEWDGqNUUNF/IrELgiA8w5ZLMVxLSGNe3wZ6P6/+NJHYBUEQCpCakcvXhyNpU8uJPo2q6DqcYhGJXRAEoQBfHY4kIyefBf0bGswUzF9EYhcEQfiXkLhHbLscw5i2NajraqvrcIpNJHZBEISnyLLMp/uv4WRtzpsveOk6nBIRiV0QBOEpZ26mcPnuQ956wQtbPT5d+jwisQuCIPxJlmW+OXqDqpUsGexbTdfhlJhI7IIgCH86GZnM1dhHTOtaBzMTw02Phhu5IAiCBv01Wq/maMmgFh66DqdURGIXBEEAjkUkEXrvMdO7emFqbNip0bCjFwRB0ABZllnyxw2qO1kxsFlVXYdTaqJsryAAqdmpJGQk8Cj7EY9yHmFhYkEV6yq4WbvhYO5gcAdUhOK5cOcBYffS+HxgI0wMfLQOIrELFZQsywTcD+B03GnOx5/nxsMbz3yti6ULHTw60MGjA22qtMHK1KoMIxXKwuozUThZmzGgHIzWQSR2oYJRqVUcjTnK2tC1RKRGYGpkSnOX5rzZ/E1q29fGwcIBO3M7svOzSchIID49nitJVzh09xC7b+7G1tSWQXUHMaz+MNys3XT9xxE04FbSE45fT+KtF7ywMDXWdTgaIRK7UGH4J/jzv4v/427aXWrY1WBB2wX0qtHrmSPwBk4NABjZYCR5qjyCkoLYdWMXG69tZOO1jbxY80VmNJ8hEryBW3M2CnMTI0a2rq7rUDRGJHah3EvLTWNxwGJ+ufkLnraeLO60mG6e3TA2KvrozNTYFL8qfvhV8SM+PZ5NEZvYfn07R6OPMs5nHGN8xmBpYqnFP4WgDSnpOewOuscrzT0MqixvYTSySiBJ0lpJkpIkSQrTxP0EQVMCEgN4+deX2XtrL+N8xrG732561OhRrKT+b+427sxpOYd9A/bR0aMjy4OX8/KvL+Of4K/ByIWy8POFaHLz1YxvX1PXoWiUppZ/1wO9NHQvQSg1WZbZErGFiUcmYm1qzZY+W3i7xdtYmFho7BlVbaqyuPNi1vZci6mxKROOTOBz/8/Jys/S2DME7cnNV7PZP5qu9Vyo42Kj63A0SiOJXZbl00CqJu4lCKWVo8rhw3Mf8vmlz2lftT1b+mz5e75cG1q6tWRn350Mrz+cLde3MPi3wUSmRmrteYJmHA5PJCU912D6mBaH4W/YFISnpOemM+XYFPbd3sfkJpNZ0nUJtmbar6dtaWLJ3FZzWd1jNRl5GQz/fTi/3vpV688VSm6zfzTVHC3pUMdZ16FoXJkldkmSJkmSFCBJUkBycnJZPVaoQFKyUhh3eBxB94NY2H4hU5pOwUgq27GLXxU/dvTdQZPKTfjo3EfMOz+PXFVumcYgFO5WUjoX76QyrFV1jIzK3+GzMtsVI8vySmAlgK+vr1xWzxUqhnvp95h4ZCIpWSl83/V7Onh00FkszpbOrOy+kmVXl7EqdBVRj6P4tvO3OFk66Swm4Z+2+Mdgaizxqq/mi309yHpAWEoYpkamWJpaYm9uT027mmV6ellsdxQMXmxaLOOPjCc9L51VPVbRpHITXYeEsZExM5rPwNvRmw/PfsiwA8P4odsP1HWoq+vQKrzsPBW7AmPp5VMFZw1tcczMy2Tf7X0cjT5KwP0A1LL6H5/3cvBicN3BvFTrJWzMtL9QK8ly6QfPkiRtBToDzsB9YJ4sy2ue9XpfX185ICCg1M8VhOi0aMYdHkeuKpeV3VdS36m+rkP6j/CUcGYcn0F6XjrfdvmWtu5tdR1ShbYrMI7ZO4PZNqk1rWuV/l1U0P0g3j/7PvfS71HTvibdq3ennXs7JEkiKy+L2Cex7L65m4jUCCxNLFncaXGJ31FKkhQoy7JvYa/TyIhdluWhmriPIBRH1OMoxh8eT746n9U9VuPt6K3rkArU0LkhW/psYcofU5h6bCqftv+Ul2q9pOuwKqzN/tHUrmyNX03HUt0nT5XHsqvLWBu2Fncbd9b2XEtLt5YFvva1eq8RnhLOjhs7ymTwIXbFCAYpOi2a8YfHo5JVrOm5Rm+T+l9crV1Z32s9zVyb8d6Z99gQvkHXIVVIN+4/4UrMI4a28izVnHeOKofpx6ezJmwNA70Gsrvf7mcm9b80dG7IgrYLcLbU/i4ckdgFgxOTFsO4w+OUpN5jDV4OhtFJ3tbMlhUvrKBH9R58HfA13wd9jyamQoWi2x0Yh4mRxMulqOKYp8pj1slZnIs/x/w285nfdj7WptYajLL0xOKpYFBin8Qy7vA48lR5rO65mjoOdXQdUrGYGZuxqOMi7PztWBW6iqz8LOa0nCPqvZeBfJWaX67co0s9lxIvmuar85lzeg6n4k7xod+HvFL3FQ1HqRkisQsGIzYtlrGHx5KtymZNjzUGu8PE2MiYj1t/jLmxOZsiNpGtyuaj1h+V+Z77iubMzRSSn+SUqp/pQv+FHIs5xpyWc3it3msajE6zRGIXDEJMWgxjD48lV5XLmh76P6deGEmSeLflu1gYW7AmbA2yLPNxm49FcteiXYFxOFqb0cXbpUTX77+zn503djLOZxwjG4zUcHSaJRK7oPf+2tKYp8rT690vxSVJEm82fxMjyYhVoaswMTLhA78PxLSMFjzKzOXotfsMb+2JmUnxf3jeeXyHTy58QnOX5kxvNl0LEWqWSOzl0eN7kJkCRiZgZAq2rmBhr+uoSiQyNZLXj76OWlazuudqg51+eRZJkpjebDr5cj7rwtZhLBkzt9Vckdw17LfgeHJV6hJNw2TlZzHr5CwsjC1Y1HERJkb6nzb1P0KhcGo1RB6AG4fg7ll4ePefn5eMwbM1ePWA+n3BqbZOwiyu4ORgJh+bjKWJJWu7r6VWpVq6DkkrJEni7eZvk6fKY1PEJixMLHi7xdu6Dqtc2RUYR/0qdjR0L/4AZ3HAYm49usWKF1bgau2qheg0TyR2Q6ZWwbVf4dRXkBwBFpWgRnvwewMqeYI6H1R5kHQNbh6BY/PgjwXQbCR0fg/squj6T/BMF+Iv8OaJN3G2dGZVj1VUtSkfTYafRZIk5rScQ44qh7Vha3Ewd2CMzxhdh1Uu3EpKJzjuMR/2Kf7BoKD7QWyP3M7IBiNpV7WdFqLTDpHYDVViKOx5A+6HQeV68MoaaDgAntUZqNvHyhTN+R/g8moI2QHt34IOs8FYv74Mdt3YxWcXP6OGfQ1Wdl9JZavKug6pTEiSxAd+Hyit/AIXY29uzwCvAboOy+Dtu3oPIwn6NXEv1nW5qlwWXFiAu7U705pO01J02qFf39FC4dRquLgM/vgELB3+TOgDwagIC0L2VaH3F+D3ujJyP/k5RJ2GQWvBVvcNmdWymu8Cv2Nd+Draubfjq05flUktdX1ibGTM5+0/50nuE+ZfmE8l80p08eyi67AMlizL7A2Op21tZ1zsitc9a03YGu48vsOybsue2fBcX4m9VYYkMxU2DYAjHyrz5ZMvQKNBRUvqT3OsCa+uhwE/QfwVWNEBos5oJeSiepj9kOnHp7MufB2veb/G0m5LK1xS/4upsSnfdv6Whk4NmXN6DmEpopVwSQXHPSb6QSb9mhZvtH7n8R1Whayid43edPToqKXotMdwE3vGA7h1TJlSuLhCmWK4/juk3AJVvq6j07zUO7CmO0Sfh75L4LVNYF3KynRNhsCEP5QdMz+/DGG7NRNrMQUkBjDot0FciL/A+37v84HfBwax80CbrEyt+KHrDzhZOjH1j6nEPonVdUgGae/Ve5iZGNHLp+jvSGVZZuHFhViYWDCn1RwtRqc9hvXdkxoF1w8ov2Ivwr9qHv/N1Boa9IOmw6B6++KPaPVNjD9sG6r8eUftg+ptNHdv1wYw8ThsHQK7xkNOOrQYrbn7P0eOKoeVIStZHboaDxsPNr24Sau9SQ2Nk6UTP77wIyN+H8GUY1P4uffPVLKopOuwDIZKLfNbcAJdvV2wszAt8nUnY0/in+jPe63eK5OCXdpgWIn97LcQtAFcfZRFv5odlblhS0cleT+4DSk3IOYihO+B4K3gUEPZAdJosGEm+BtHYMdIsHOH4bu0s1XRwk65946R8NsMyE2HNlM1/5ynnIk7w+eXPif2SSz9avfjfb/39a6Qkj6oaV+T77t+z8QjE3n75Nus7L4SU+OiJ6mK7PztFFLSc+hfjGmYPFUeiwMXU9O+Jq96v6rF6LRMluUy/9WiRQu5RB7ckeXUqKK9NidDloN3yPJPnWR5np0sr+goy1FnS/ZcXQn/VZYXOMnyig6ynJ6s/efl5cjy9pHK39fFFVp5xNWkq/LUY1Nln/U+ct89feWL8Re18pzyZt+tfbLPeh95/vn5slqt1nU4BmHWjquyz8eH5Kzc/CJf83P4z7LPeh/5VOwpLUZWckCAXIQca1gjdseaRX+tmRU0fhV8XoGwXXBsPqx/EZqPhp4LwVz77alKJXgb/DoZPFrC8J0lPjmqltU8zH6o/Mp5SFZ+FkaSEUYYYWpsioO5Aw4WDtib22NiYgavrAX1aDg4B8ztoGnpe6jkqnI5d+8cG69tJOB+APbm9rzV/C1GNRglRp9F1Ld2X24/us2asDXUqVSH4fWH6zokvZadp+JwWCK9fNywMH3GFuB/eZzzmB+Df6RNlTZ0qKq7nrmaYFiJvSSMjKDxYKj3Epz6As59D3fPwMBV4FFohyndCNoI+2ZAzQ4wZGuxfgilZKVw9t5ZQpJDiHwYyc2HN8nKzyr0OmPJGFcrV9xt3HF3r4ZnTmOqH51Ntfw0PBoNwc7MrsjH3GVZJiEjgdCUUE7GnuRk7EnS89JxtXLl3ZbvMtBroMFtH9MHM5rP4M7jOyy6vIiadjVpW1W02HuWUzeSeZKTT99i7F1fEbyC9Lx0ZrecbfAlHTTS87S4dNrz9O5Z5WBPWrxyaKfdm6BP/4gBa2H/21C7GwzZDKaWhV6SlJnE3lt7+SPmD8IfhANKUwdvB2+8Hb3xtPXE0cIRBwsHrEysUKNGlmWy8rN4nPOYhzkPSc5MJiEjgfj0eOLS40jKTPrHM2xMbXC3ccfZ0vnve5kZmSmjf8mItNw0HmU/4kH2A248vMGjnEcA2JnZ0dWzK92rd6dNlTZihF5KmXmZjDg4gvsZ99n20jaq2VbTdUh6acbWK5y5mczlD17AxLjwtbX49Hj67OlDv9r9WNB2QRlEWDJl2vPUoNRoD2+cVRYJj82DewHQf7mygKhrl1bB77PBqycM3gimzz5QIcsyF+IvsD1yO6fiTqGSVTSp3ITpzabT0aMj3g7epRp1ZOVnEZscRuzeN4jLe8y9Rr2IV2WRmp1KdFo0qdmp5KvzUckq1LIaW1NbHCyUaZ0u1brQwKkBDZwaUN+pPqZGIplripWpFUu6LGHI/iG8deItfu79s3j38y/ZeSqORdynf9OqRUrqAMuvLscIIyY3mazl6MpGxUvsAJaV4NUNcGEpHJ0Hyd2UfeGVdVgO9vwPysEj7z7w6jowKbjDiyzLnL13lh+DfyQ0JRRHC0dGNxzNK16v4GnnqbFwLE0sqVulJXWH7IHVL0DQ7zDhWIH1ZWRZLtO3rrIsE/cwi9B7j4lKycDC1Bgbc2MqWZnhV9ORSlZmZRaLLlSzrcaijouYfGwy88/P58uOXxr81IEmnYxMIjNXxUuNi1YL6c6jO/x25zeG1x+Om7XuT2BrQsVM7KBMv7SdDlWawq6xsKobDFgB9cu4e7wsw4nP4PRXSq2XASvBpODEFJocypeXvyQ4OZgq1lWY12Ye/Wv31+70RiVPZfF23YuweRCMPfifdzdllVTupmSw5mwU+0PieZiZV+BrjCRo7ulA1/ouDGnpiaN1+Uzy7aq2Y0bzGSwJWkJD54aMblg2Zw8Mwf6QBJyslR/yRbH06lIsjC2Y0GiCliMrOxU3sf+lZgeYdBK2j4Ttw6HjHGXfe1nseVer4dBcuPQTNB8FL31XYBGvB1kPWBK0hD239uBs6cxHrT9iQJ0BZTdfXaUJDN4AmwfDjlFKoi/DufLIxCd8e/QGh68lYmpkRO9GbrSs4UhjD3vquNiQm68mI1dF4uMsTkUmczwyiUWHIll2/Bbj29dkfIda2FuWv+mg8T7jCU8J59vAb/Fx9qGFawtdh6RzWbkq/ohIYmDzok3DhD8I52j0Ud5o8gaOFkX7QWAIKt7i6bPkZcOBWXB1k7JwOXAlWGvx1FluhrKIG7EP2kyDHv/7zyKuLMvsvb2XRZcXkZWXxYgGI3i98evYmOloq2bQz7BvGjQdAf2Xan3RWaWWWXn6Dt8evYGFqREj21RndJsaRSrmdPP+E747dpMDoQnYWZjwYZ8GvOrrUe6mLNJz0xlyYAgZeRns7LvTYE9KasrvoQlM2RzElol+tK1d+N/FG8feICwljIMDDxpEbaKiLp4a4FFMLTG1UJLVS98qO2dWtFfqsmjD43uwthdc36/sqS8gqSdmJDLljyl8dO4jvCp5sbv/bmb5ztJdUgdoPhI6vav88Du1SKuPin6QwasrzvPloet0q+/CidmdeadnvSJX6PNytWXZ8OYcmNGe+lXsmLM7hDc2BZKakavVuMuajZkN33T+hvTcdN459Q756nJYJ6kYDoQk4Gxjhl/NwusoXU26yrl75xjnM84gknpxaCSxS5LUS5KkSEmSbkmSNFcT99QJSQLfccoioaklrH8JTiyE/BzNPSPqNKzqotS9GbZDObr/r6R+MOogA/cOJPB+IHNbzWVdr3XUsteT7kGd34MmQ+HkQmUErwUBd1Ppv+wct5LSWTKkKcuHN8fJpuDF5MI0dLdn68TWvP9iPY5fT6Lnd6fxv/NAwxHrVl2Hunzc5mMC7gfww5UfdB2OzmTm5vPH9fv09qmCsVHh78yWXV2Go4UjQ7yHlEF0ZavUiV2SJGNgGdAbaAAMlSTJsCs5VWkMk04pJXFPffnn6P1C6e6Zkw77Z8KGvmBmAxOOglf3f7wkMy+TD89+yJzTc6hVqRa7+u5ieP3h+tW5XpKg7/dQu6uyZfT67xq9/cHQBIat9sfByozfprenf9OqpZ4+MTKSmNSxNnuntsfWwoQRa/zZFRinoYj1Q9/afXm17qusDVvLqdhTug5HJ05cTyY7T03vRoXvbAm8H8jFhIuM8xlXLreLaiJjtAJuybJ8R5blXGAb0F8D99UtCztlnn34LmX+fV0v+GUSJEUU7z6qfKUc7o9tlMNHbaYp++hd/tmm63rqdQbvH8y+2/uY1HgS63ut1+j2RY0yMYPBP///jqKYixq57cYLd5myJQgfdzt2T25LdSfNFgVr4G7HnsntaFnDkdk7g/n6cCRqddmvMWnLu63epZ5jPd4/+z7x6fG6DqfMHQz7azdM4dMwy68ux9nSmcHeg8sgsrKnicReFXi6WHTcn79XPnh1h6kXlROqEb/B8tawZQjcPKosgD5LZqpy4OiH5rBrHJhYwrhD0PMzpY7NU/bc3MOI30eQlZfFmp5rmN5suv7XIze3UXbH2HvAlsFwP7xUt/v5YjQf7w3nhfqubJnYWmvbFO2tTNkwrhVDWlZj6YlbzN4ZTL7qGeWfDYy5sTmLOy1GJat459Q75KkK3hJaHmXnqThxPYkeDV0LnYa5lHCJS4mXGO8zHkuTwk92GyJNZI+C/hb/MwySJGkSMAnA01NPR6LPYmYN3T+Bdm/BpZXgvwJuHAQjE6VIl0t9MDZXtgBmpEDcJXhwS7nWo5WSzL1f/M9Wxuz8bBb6L2TPrT34VfHjyw5f4mRZyuYZZcnaGUb8oiwEb+gHYw6AS71i32bH5Vg++jWMbvVcWDasOWYm2p16MjU24vOBjahayZLFR2+QladiyZBmWn9uWfC08+STtp8w69Qsvgn8hndbvavrkMrEmZspZOSq6O3z/ENJsiyz7OoyXCxdGFR3UBlFV/Y0kdjjgKcLVngA/3kfKMvySmAlKNsdNfDcsmflCJ3nKqP3mAtw55SyGBr+K6jzQZUH5rZKsm86HGp2Ao+C9xYnpCfw1sm3uPbgGhMbTWRq06kYP6sRtT5zqA6jf4P1fZT1g7G/g7NXkS/fe/Ue7/4SQgcvZ5YN135S/4skSUzv5oWVuQmf7r9G9s8B/DiiRZErAeqzHjV6MOz+MDZFbMLX1Zdu1bvpOiStOxiWgL2lKW1qP39g5J/oT1BSEO+1eg8Lk+L1QDUkpd7HLkmSCXAD6AbcAy4Dw2RZfuZ7c73cx16GLideZtbJWeSp8/i8w+d0rtZZ1yGVXvINpSyykQmM3g/OdQq95MzNZMauu4xvDQfWjWmFpZlukuoW/xg++DWU9nWcWTXKt1wk91xVLqMPjiY6LZodfXfgYeuh65C0Jjdfje//jtK9gRuLBzd55utkWWbUwVEkZCRwYOABzI1LttNKl8psH7ssy/nANOAwEAHseF5Sr8hkWWZzxGYmHplIJYtKbOmzpXwkdYDKdZW2fao8WNsT4q8+9+XX4tOYvCmIOi42rBzlq7OkDjDMz5MvX2nMmZspTN4USE6+SmexaIqZsRlfdfoKJJh9aja5qvK1f/9p52+nkJadT+9C+ppeiL/A1eSrTGw00SCTenFo5H2vLMu/y7JcV5bl2rIsf6aJe5Y3uapc5l+YzxeXvqCjR0e2vLiFmvbFaBxiCFwbwLjDYGqlnAGIOl3gy+49ymLs+kvYWpiwbmzLYvWj1JbBvtVYOKARJyKTmbr5Crn5hr+g6mHrwaftPiX8QTiLAxbrOhytORSWiLWZMe29nn3SVJZllgUvw83ajQFeA8owOt0w/NUiA5CSlcL4w+P55eYvTGo8ie+6fKfbE6Ta5FwHxh9WdstsegVCd/3j00+y8xi37jKZuSrWj21FFXv92ZUwzM+TT/o35FjEfd7afqVc7Jbp5tmNEfVHsOX6Fo5GH9V1OBqXr1Jz5Np9utV3fe4U2l/NZyY1noSZcfksDPc0Pd9TZ/giHkQw48QMHmU/4qtOX9GrRi9dh6R9du7KIuq24bB7PCQEQ7d5qCRjZmy9wq3kdDaOa4W3m/4d4x7Vpga5+Wr+dyACS9NQvhrUGKMinGLUZzNbzCQkOYSPz32Mt4O3ds9H5OfA7eOQHAlPEpSGNpKRsvHAyklpLu/eDJy9wbj06efS3VRSM3KfOw0jyzI/Bv+Iu7U7L9d+udTPNAQisWvR0eijfHD2A+zN7dnYeyP1neoXflF5YeUIo/bC4ffg/PeQGMI3tnM5EZnKZwN8aFdHf4tVTehQi8xcFd8cvYG1uTEL+jU06OJhpsamfN3pa17d/yqzTs1i04ubND/HHHsJrvwM4Xsh57Hye2a2/1+/P/MBZD0E+c93QSaWSmvKOi8oZ0VcGpSoqNyR8PtYmBrRybvyM19zMvYkoSmhzG8zv8J08BKJXQtkWWZFyAqWX11O48qNWdJlScWsumdiBn0WQ5WmqH57m5Hqobj5fMRwvz66jqxQ07vWISMnn59O38HKzIR3e5WuI5WuVbGpwsL2C5n6x1S+uPQF89rM08yNs9PgyAdKn15Ta6jfV2kiX81P2fr7NLVKqZEUfwXigyDqjNLF7Ng8sPcEnwFK83m3xkVK8mq1zKGwRDp6VcbKrOBUppbV/HD1B6rbVadfnX6a+BMbBJHYNSwrP4uPzn3E4buH6VurL/Paziv3K/CFOW/fmy9y0/jR6idG3poJByKg+wLl4JeekiSJub3rkZGbz4pTt7G1MGFql8K3cOqzjh4dGe8znjVha2jm0ox+tUuZ6O6cgr1TIe2ecraj07vP/zc1MlbWYJzrKMkflKmaW8eUU90XlsG5JeBcF5qNUIrN2bg883Yh9x6TmJbNHJ9ndz47FHWImw9v8mWHLytUi0aR2DXofsZ93jzxJtejIVoOAAAgAElEQVQeXOPtFm8ztuFYgx7laUL0gwymbA7C2ckH20nn4OzncHEZ3DwCvReBt/6uOUiSxCf9fMjMUfHV4UiszIwZ286wdzJNazaNkJQQPr3w6d/NzkskfI9SKsOhprITqlqrkt3Hzl1pMtN8lFKG49peCN4GRz+GPz6Bur3A73Wo0eE/o/hDYYmYGEl0q+da4K3z1Hksu7oMLwcvetXU368zbRC7YjQkLCWMoQeGEvU4iiVdljDOZ1yFT+pPsvMYv0E5iLZ6lC92NrbQa6FSesDEAra+BluHwcNoHUf6bEZGEosGNaZnQ1cW/HaNHZdjC79Ij5kYmbCo4yJszWyZeXImT3KfFP8mN4/B7olKuYzXT5c8qf+blSP4jlV2VU29DK2nKCe8N/SFH9tC4HqlIB/KdOfh8ETa1HbC3qrgkfi+W/uIeRLD9KbT9atCahmoWH9aLTkYdZAxh8ZgamTKzy/+TBfPLroOSedUapkZW69wNyWD5cObU8P5qbfoNdorFS5fWAB3TsDSlnBsvjJfq4dMjI34fmgzOng58+4vIey9ek/XIZWKs6UzizsvJj49ng/PfohaLsa2zugLsH2EUhNo2HalGJw2VK4LPT6Ft8Oh/zKQjOG3N2FJEzi/lFv37hOVkkGvZ+yGyVHlsCJkBY2dG5efQ4DFIBJ7KahlNT9c+YE5p+fQ0KkhW1/aSl2HuroOSy98diCCE5HJzO/XsOAWZSZm0P4tmHZZaeJ99lv4vhlcXqOUOtYz5ibGrBzpS8sajszcEcyhsERdh1QqzVyaMdN3Jsdjj7M6dHXRLkqNgi2vgX1VGLEHLCtpN0hQGt40GwFvnIGRvyoJ/8gHeKxvxSST/XT3sivwsi0RW0jMSGRG8xkV8p2zSOwllJmXycyTM1kZspIBdQawusfqctUMtzQ2XYxm7bkoxrarwYjW1Z//YnsPGPgTTDwBlb3hwExY0U55u69nLM2MWTumJY2q2jN9axAnI5N0HVKpjKg/gj61+rD0ylJOxxV8SvhvajXsmw7ISkVPm2dvL9QKSYLaXZSCc+OPEU5t3jfZgsv6dsqOHPX/l4F4mP2QVSGr6OjREb8qfmUbp54Qib0EEtITGHVwFCdiTzCn5RwWtF1QYfbHFubszRTm7Quni3dlPuxTjEZaVZsrc++vbQZVLmx+RTm5mnpHe8GWgI25CRvGtsLLxZbXfw7k3K0UXYdUYpIkMa/NPLwdvZl7ei7Rac9Z6whcB3fPKP15HQr5Ya1lsdYNGZQ+mwPNV4FtFeUHzspOf3c5+ynkJzLzM5nVYpZO49QlkdiLKeh+EEMODCE+PZ7l3ZYzssHICvlWryA37z9hyuZA6lS24fuhzYrUd/IfJAnqvwRT/JUm3zH+sLwNnP4a8vWniJW9lSmbJvhR09ma8Rsuc+G24fZQtTSx5Lsu32FsZMybx98kPTf9vy96FKPsUqnVWdm9omOHw5VpsMbtX1L6Ew9ap+yoWdeLuztHsP36Nl7xeoValfSkT7AOiMReDLtu7GL8kfHYmdmxuc9m2lVtp+uQ9EZSWjZj1l3G3NSY1aN9sS1NYS8TM6XJ97RLULcnHP9UGZEVty2hFjlam7Fpgh/VHKwYv+Eyl++m6jqkEqtqU5WvO33N3bS7zDk9B9VT0xrIsrJoKctKr1s9GMQcCkukobsd1RytlHh8BiprNR1m812KP2aqfCabuCkxV1AisRdBnjqPhf4LWXBhAX5ufmzus7n8VWYshYycfMZtuExqRi5rR7dUvuE0wc4dBm+EodshIxlWdoErm/TmG9bZxpzNE/1ws7dgzNpLBBhwcver4sf7fu9z5t4ZFgc+VQkyYp9S+6X7Ap1PwYAygAiMeUjPhv/aDWNmzaUGPfnDyoLxsi3OB2bDpoHwyLC3p5aUSOyFeJj9kNePvs7W61sZ1WAUS7stxc6s4JX4iihfpWb61itci09j2fBmNPKw1/xDvHsp2yOrtVROOu55Qyk2pQdcbC3YOrE1rnYWjF57icBow03ug70HM7z+cH6+9jO7buxSFkxPfglOdcB3nK7DA+DItfvIMv/Z5pinyuN//v/Dw8aDUSP+gN5fKVN5P7aFq1v0ZjBQVkRif47I1EiGHhhKcFIwC9sv5J2W7+h/k+kyJMsy7+8J5fj1JD7p70PXZ5wA1AhbN2W7W+f3IGSb0kA7p4D5YB1wtbNg66TWuNhZMHrtZQKjH+o6pBKb7TubdlXb8dnFzzjv/y0khSulAvSkbePh8ERqOVvj5fLP/fPrw9cT9TiKD1p/gIWZFfhNgsnnwNUHfp2s7L3PMNyF7uISif0ZDt89zMiDI8lT57Gh9wb61u6r65D0zqLDkewIiOPNbl6Fb2vUBCNjpefsyz8qBaQ29lcWzfSAq50ycq9sa/7nyN0wk7uJkQlfd/ya2pVq83bkeq5VrqUU5tIDjzJzuXD7AT193P6xYSHuSRw/hfxE9+rdaV+1/f9f4FgTxuyH7p8qJSyWt4Gb5a8mfUFEYv8XlVrFd4HfMfvUbLwdvNn+0nZ8nH10HZbeWX3mDj+evM1wP0/eeqHozas1oukwZe49MQTWvag3yd3NXknuzjZmBp3cbcxsWO75MvaqfKZUMic24z+96XXij4gk8tUyvZ6aX5dlmc8vfY6xZMyclnP+e5GRMbSboZyTsHaGzYPg9zmQl1WGkZc9kdif8jjnMdOOT2NN2Bpe8XqFNT3XVMxyu4XYFRjH/w5E8GIjNz7p76Ob7Z71X4LhOyH1NmwfqTfbId3sLdg2qY1hJ3e1GpcLK1iRbUmekTGTj03mQZbut3QeCk+kir0FjZ9axzl09xCn404zpekU3Kyf0/PUzUdJ7n6T4dJPsLKz0gCmnBKJ/U83H95k6IGhXEy4yEetP2J+2/kVooVWcR0ISWDOrmA6eDnz7WtNi79XXZNqdVbqiESfhf1v6c0CmcEn95tH4H4otdq/y9JuS7mfcZ9JRyfx+K8GGjqQmZvP6RvJ9Gz4/9MwCekJfHrhUxpXbszw+sMLv4mpBfT+Qjk5m/0YVnVVzkioDb95+b+JxI7S6Wj478PJys9ibc+1DPYerOuQ9NKJ60m8ue0KLao78NPIFpib6MGCWuPB0GkuXN2s1JvRE/9O7kExBpTcA9aCjRv4DKSZSzOWdFlC1OMoXj/6esmqQWrAychkcvLVf++GUctqPjj3AflyPl+0/6J4mxrqdIPJ55WmIMc/hTXd4X64liLXjQqd2FVqFUuCljDz5Ey8HLzY/tJ2mrk003VYeun87RTe2BRI/Sp2rBnT8pkda3Si81zwGQR/LFD2XOsJN3tlt4yTjRmj11ziauwjXYdUuEexcOsoNB8Jf5bJaFu1Ld90/obI1EimHJtCZl5mmYd1MCwRJ2szfKs7ALAxfCOXEy/zXqv3qGZXrfg3tHJUTqy+sgYe3oWfOsLxz7S/jVZdNg3SK2xif5zzmKl/TGV16GoG1R3Eup7rcLF6dreWiuzy3VTGrw+gupMVG8a1wq40p0q1QZKg/1Kl887eacrbbD1Rxd6SrRNb42Btxsg1/oTG6U9sBbryszKl9a/SAZ2rdWZRp0WEpoQy8cjEMp2Wyc5TcTziPj0aumJibET4g3CWXFlCN89uvFynFM2pJQkaDVJqv/sMgtOLYFkrpYmIpqf18nOUd0I/NIOk65q9dwEqZGKPTI3ktf2vcSnxEvPazGNem3liPv0ZrsQ8ZOy6y1SpZMHmCa1xtNbTvydTS3h5BTxJhEPv6Tqaf3CvZMnWSa2xszBl5Fp/rifqZ915VPkQ9LPSYLqS538+3b16dxZ3XkxEagRjDo0hOTO5TMI6czOFjFwVvX2qkJiRyPQ/plPZsjLz2szTzMK9tZNSYXTkHqVv684xsPoFuHOy9Ak+Ow0u/ghLmsL+t8HKCcrgHU+FS+wH7hxgxO8jyFPlsa7XOgbVHaTrkPRWaNxjRq29hJONGVsmKHu09ZpHC2j/tjLffv13XUfzD1UrWbJloh/mJkaMWO3P7WT9OFz1DzePwJN4aDHmmS/p5tmN5S8s5176PUYfGk1smvaP7B8MTcDe0pQmnpZMPz6dzPxMlnZbioOFg2YfVLurUve9/zKlj+vG/rC8NVxaVbwmMGq1cur116mw2BsOzVX21I/cAxP+UCqZapkkl+InkiRJrwLzgfpAK1mWA4pyna+vrxwQUKSXakyeOo9vAr5hU8QmWri24OtOX4utjM9xLT6NYasvYm1mwo432lC1kqWuQyqa/FxY1QXSk2CqvzKXqkduJaUzZOUFTIyM2PlGG83V1dGEzYOVLYBvh4Px89dQgpODmXJsCpIk8U2nb2hVRUPt8f4lN19Ni/8dpUcDF3Kd1nL63mmWdVv2z4NI2pCXDeG/wKWVEH8FjEygqi/U6qR8tKkM1pWVFo/pSZCRBA9uKQfn7p6BzAfK6L/RK9B8NHj4aiQsSZICZVku9GalTez1ATXwEzBbXxN7SlYK75x6h4D7AYyoP4KZvjMrVMfy4rp5/wlDVl7EzMSI7ZPa4OmkR8mnKBJDlcWwVq8r29v0zPXENF776SIOVqbsfKOtfrwTehQLSxpDh1nQ9cMiXRKTFsP049OJTotmbqu5DKk3RONhnYhMYuy6i3Ruf5zAB3/wvt/7DK03VOPPea64QLj+G9w5BQlX4XmtBO2qQs2Oylbcen3A3FajoRQ1sZdqa4MsyxF/Pqw0t9Gq0ORQ3jr5Fmk5aSxsv1CUBihEVEoGw1b7Y2QksWVia8NL6gBujZTFv8uroNVEcKqt64j+oZ6bHWvHtGTEan/GrLv09/y7ToVsVxJWs5FFvsTTzpPNL27m3TPv8pn/ZwQnB/Oe33saLZL3W0gUNtV/JvDBdaY3m172SR2UKT6PFsp/Zz2C5EjITFEqjubnKCN3G1elZWCl6npR2rhUI/a/byJJJ9HDEfvuG7v5zP8zXKxc+K7Ld9RzrKf1ZxqyuIeZDF5xgZx8NdsmtcbLVbOjjTL15L7SQ7VON3jtZ11HU6CTkUlM2BBAi+oObBjXCgtTHZ4L+LE9mFnD+MPFvlSlVrEyZCU/hfyEs6Uzn7T7hLbubUsdUlLGA7ptGg0WMcxr87FYD6PoI/ZCF08lSTomSVJYAb/6FzOgSZIkBUiSFJCcXLLV9KCYhxwOTyQ77/knxXJVucw/P5/5F+bj6+rLtj7bRFIvRFJaNiNW+5Oek8/G8a0MO6kD2LoqzbIj9kHMRV1HU6DO3i4sHtwE/6hUZu0MRq3W0cnZlJtwP1RpKl4CxkbGTG46mc0vbsba1JrXj77O+2feJ+5JXInuJ8syB6MO0v/Xl5HN7jGqlkjqxWVQI/aZO67yS9A9LE2N6VKvMn0bu9Ojods/jrUnZiQy8+RMQlNCmdBoAtOaTsNYT0qO6quHGbkMWXmR2IeZbJrgR3NPDe820JXcDPihhdIwe/xRvXiLXJCfTt3m84PXeaNTbeb21sEA5NQiOLEQZl5TmpuUQnZ+NiuCV7ApYhMqWcUrXq8wpuEYPGw9inT9nUd3+CbwG07FncLeqBap0f0JnDMaSzPxPQxlNMde1r58pTEDm3lwMCyBw+H3+T00kZrO1kzuVJuXm1UlOCWQ2admk52fzbedv+WF6i/oOmS9l5mbz9j1l4l6kMH6MS3LT1IHZWqhywewbxpc368cIddDkzrWIiY1kxWnblPN0ZLhfmXcqSh8D3i2KXVSB7AwseCtFm8xrP4wVoasZPeN3WyP3E4j50b0rNGTVm6tqGpb9e95+My8TBIzE7mUcIl9t/cRmhKKpYkls1q8w7e/uPCCl6tI6iVQ2l0xA4AfgMrAI+CqLMs9C7tOE3PsKrXM0WuJLD1xi7B7j6lc1Z9cu3142nmypMuSCt3ItqjyVGombAjgzM1kVoxoQY9/txsrD9QqZdRu5ajsIdbTUXu+Ss3EjQGcvpnC2jEt6VS3ctk8OOk6LPdTOg75TdL47RPSEzh49yCHog4Rkfr/PWttTG0wNjL+xwnWug516Ve7H31q9SEsRs2YdZdZObKcfl2WUFntitkD7CnNPUrK2Eiil08VOtS1Y8qRuQQ9OEleWkOqm82gkmnR3vZVZGq1zDs7gzl1I5kvBjYqv988RsbQdjocmAnR56CGlvc/l5CJsRFLhzXnlR/PM31LEHuntaems7X2Hxy+B5CgQbGWzIqsik0VxvmMY5zPOGLTYol8GMm99HvEp8ejklW4WbvhZu2GVyUvvB29/75uf0gwtuYmdPIuox9w5YxBTcX8253Hd3j7xNvcTbvLtKZvknG/Az+eus2FW6f4YWgz2tYRB5Ce5YtD1/n1ajyze9RlSKv/Hh8vV5oOg5Ofw9nv9DaxA1ibm7BqlC/9l51jwobL7JnaTrvbIGVZSew12iuLzVpWza5akQp25eSrOByeSI+GbvpRQdQAGWxJgcN3DzN0/1Ae5TxiZfeVvN5kAjN7eLN/egccrc0YufYSa85GoYnF4fJm08VoVp6+w8jW1ZnapY6uw9E+U0vwe12pWpgYputonquaoxXLhzcn+kEmb227ikqbO2WSIiAlEhqWopCWFpy+kcKT7HxealJF16EYLINL7LmqXD67+BmzT82mjkMdtr+0Hb8qfn9/3tvNlj1T29Gtnguf7r/G7J0h5OSXv0L6JXUiMomP94bRxbsy8/o20OvDZRrVcoJyxPv897qOpFCtazkxr19Djl9P4tujN7T3oIh9gAT1+2nvGSWwPySeSlamtBfvuEvMoBJ7bFosI34fwbbIbYxqMIr1PdcX2A7LxtyEFSNa8GY3L3YHxTFpYyBZuSK5RySkMW1zEPXc7Fg6rDkmxgb1z186lg5KcavQXfAoRtfRFGpk6+q85luNpSducfz6fe085MYhqNYKbPSnXHVWroqj1+7T28cN04r09alhBvU3tyJkBXHpcSzpsoR3Wr6DqfGz5x+NjCTe7l6XL19pxOmbyYxZd4n0nPwyjFa/pKTnMGFDALYWpqwd0xJrc4NeXimZNlOUj5dX6zaOIlrQvyENqtjx9vZgYlM1XOr1SaJS3KpuoZvYytSJyCQyc1X0bVz6rZcVmUEl9rmt5rKz7066enYt8jWvtfTku9eaEhD9kOGr/XmclafFCPVTbr6aKZuCSEnPYdUoX9zsLXQdkm7Ye4B3b7i6RW+aXz+PhakxK0a0QC3LTN4cWOiJ62K5eUT5WLeX5u6pAXuv3sPZxhy/Wk66DsWgGVRitzWzpapN1WJf179pVX4c3pxr8Y8Zv/5yhZqWkWWZefvCuXQ3lUWDGtPoqQ7vFVKLsUrxpsgDuo6kSDydrPhmcFPC7qXx6f5rmrvxjcNg5wEuDTR3z1J6mJHL8etJ9G/qrtsm6eWAQSX20ujR0I0lQ5oRFPOQyZsDyc0vm96DurbpYjRbL8UwuXNt+jct/g/Fcqd2F7D3hMD1uo6kyLo3cGVSx1ps9o/hYGhC6W+Ylw23TyjTMHq0eL4/NIE8lczA5uLrtLQqTGIHeLFRFRYOaMTJyGRm7QzW7lYyPRBwN5UFv12jaz0XZvfwLvyCisDIWCnpe+ckpN7RdTRFNruHN0087Hl3dwhxD0s53x59FvIy9G4a5pegOLxdbWlQRXNlfyuqCpXYAYa08uS93vX4LTiezw5EFH6BgUp6ks2UzUF4OFjy7WtNxVvbpzUbDpIxBG7QdSRFZmZixA9DmyPLMGPrFfJUpXjHeeMwmFhCzQ6aC7CUolIyuBLziIHNq1acLbhaVOESO8DrnWozrl1N1p6LYuOFu7oOR+PyVGqmbb7Ck+x8Voxsgb2l6Bb1D3buymj16maDWET9i6eTFQsHNiIo5hHfHSvh/nZZVrY51uqsHNzSE3uC4pAkxHShhlTIxA7wQZ/6vFDflfn7wrW3T1hHPv/9OpfupvLFK42o5ybe1haoxRiDWkT9S98m7rzmW43lJ2/jf+dB8W+QHKns49ejbY5qtcwvV+7Rvo5zxd2xpWEVNrEbG0ksGdKUBu52TNtyhfD4x4VfZAD2h8Sz9lwUY9rWEKOf56nTDWzdIXi7riMpto/7NqC6oxUzdwQXf/vujUPKR68emg+shAKiHxL3MIsBzcTXq6ZU2MQOStGlNaNbYmdhysQNASQ/ydF1SKVyK+kJ7+4KoUV1B95/sb6uw9FvRsbgM1CpH5OZqutoisXa3IRvX2tKYlo2H+8tZu2b23+AS0OlP6ee2B0Yh5WZMT3La4VRHajQiR3A1c6CVaN8Sc3M5Y1NgQZbVyYjJ583NgVhYWrMsmHNMTOp8P+0hWs8GNT5cO1XXUdSbM08HXizmxd7r8bz65V7RbsoN1NpE1i7i3aDK4a07Dz2BcfTt7F7xTwNrSXiux9o5GHP1682ITD6IR/uCTO4ipCyLDP3l1DuJKfz/dBmYp6yqNwag3NdCNmp60hKZErn2vhWd+CjvWHce5RV+AXR50GVq1eJ/dcr98jKUzG8dTkvHV3GRGL/00uN3ZnRtQ47A+NYczZK1+EUy/rzd/ktOJ5ZPbxpJyriFZ0kQaPBEHPeIAqD/ZuJsRHfDG6KWi0ze0cRmmHfOQHG5uDZtmwCLIQsy2zxj8Gnqh2NPSrpOpxyRST2p7z1Ql16NXRj4e8RnLqRrOtwiiTgbiqfHYjghfquTO5UW9fhGJ5Gg5SPYbt1G0cJeTpZMa9vQy7cecDac4UMSG4fB8/WYGZVNsEVIijmIdcTn5R9j9cKQCT2pxgZSXzzWhO83eyYtiWI28npug7puZ4+hLR4cBOMxCGk4nOsCR4tDXY6BuBVXw+6N3Bl0aFIIhOfFPyiJ4mQdA1qF72AnrZt9o/BxtyEfk1EJUdNE4n9X6zMTFg1qgVmxkZM2BDA40z9rAaZp1IzbcsV0rLz+HGEOIRUKo0GQ1I43A/XdSQlIkkSnw9shJ2lCW9uu1LwBoDbJ5SPepLYH2Xmsj8kgQHNqopFUy0Qib0AHg5W/DSyBXEPM5myJbB0x7e15JPfrnEpKpUvBjamvqitUToNByglBkJ36TqSEnO2MWfRoMZcT3zCN0cKOJV65wRYOYOrT9kHV4BdgXHk5qsZ5icWTbVBJPZn8K3hyBcDG3Pu1gO92ymz2T+any9G83rHWrwsDnWUnk1lpW7Ktb3KkXsD1bWeK8P8PFl55g4Xnz6VKsvKiL12FzDS/bd8vkrNxgvRNPesJAYlWqL7f2U99koLD2Z082J7QCw/nrqt63AA8L/zgHl7w+nsXZk5verpOpzyo34/SL2tzEMbsA/71KeGkzWzdgSTlv3nNOL9cMhIglr6sc3xYFgiMamZTOooFvu1RST2Qrz9ghf9m7qz6FAkvwXH6zSWqJQMJm8OwtPJiiVDmomKjZpU7yVAgmv7dB1JqViZPXUq9dc/T6XePq581IP967Iss+LUbWpVtqZHA1ddh1NuicReCEmSWDSoMa1qODJzx1VORibpJI7kJzmMWusPwJrRLcViqabZuoJnG4gw7MQO0LRaJWZ09eLXq/HsuRKn1J539laqWurY2VsphMen8XrHWmIXlxaVKrFLkvSVJEnXJUkKkSRpjyRJ5fKUgbmJMavH+FLX1ZY3NgVy+W7Z1hZJz8ln7PpLpDzJZc1oX2o6W5fp8yuMBv2UqZiUW7qOpNSmda1DqxqOLNhzFXX0eajVSdchAbDi1G1c7czF2pCWlXbEfhTwkWW5MXADeK/0IeknOwtTNoxrhbu9JePWXSbsXtlUg8zOUzF5UyARCU9YPrw5zTwdyuS5FVL9vsrHcjBqNzaS+HZIU5oY3cYoP4s8T9031QiNe8y5Ww8Y164m5ibGug6nXCtVYpdl+Ygsy/l//u9FwKP0IekvZxtzNk3ww87SlGGrLhIYrd2Re1auiokbAzhzM4XPBzaiSz0XrT6vwrP3gKotykViB6hayZIP6yehliWWRum+cuLyk7ewNTcRWxzLgCbn2McBBzV4P73kXsmSbZNa42htxvDV/lqbc0/PyWf0ukucu5XCV4MaM9i3mlaeI/xL/X4Qf8Uga8cUxCvzCvFWdVlyLlmnDWWCYh5yMCyRse1rYmsh1oe0rdDELknSMUmSwgr41f+p13wA5AObn3OfSZIkBUiSFJCcbBh1WJ6lmqMVO99oSy1nGyZsCFAWqDQo6Uk2w1f7ExT9kCVDmvGqSOplp0E/5WPEb7qNQxNyMyH2Em5NetCgih1vbw8mNrWUjbBLQJZl/rf/GpVtzXm9Y60yf35FVGhil2X5BVmWfQr4tRdAkqTRwEvAcPk5p3hkWV4py7KvLMu+lStX1tyfQEcq25qz7fXWtKjuwNvbg3l/TyjZeaWv5X7h9gP6fH+WyMQ0fhzRgr6ijkbZcqylnM6M2K/rSEov9iKo8zCp05kVI1qglmUmbw7UyNdpcfwemkhQzCNmda8rygeUkdLuiukFvAv0k2W57IcCOmZnYcqmCX680ak2W/xj6Lf07LOLMBUiX6Vm2YlbDF99EVsLE/ZObU93sc9XN7xfVJJiRgl6iuqTO6fAyAQ82+DpZMU3g5sSdi+NBb+VXU2cnHwVXxyKoJ6brXjnWYZKO8e+FLAFjkqSdFWSpBUaiMmgmBobMbd3PTaOa0VqRi4vfn+GObuK/pZXrZbZHxJPj29P89XhSF5q7M6+ae3xdrPVcuTCM3n3BlkNN4/oOpLSiTqtVK40U7bHdm/gypTOtdl6KZZ1hZX41ZCN56OJTc3i/RfriwN1ZahU74tkWa6jqUAMXce6lTn0VkeWnbjFZv8Y9ly5R/+mVeni7UKb2k44Wpv9/VpZlolIeMLZW8n8eiWeawlp1HW14aeRLejRwBVJEt8AOuXeDGyrQOQBaDpU19GUTNYjSLgKHd/5x2/P7uHN7eR0Ptl/DQ8HK62+K4x+kMF3x27Q2bsyHesa/vSrIRETXhrkbGPOvL4NmdSxFstO3OLXK/HsCoxDkqxyi3kAAAyqSURBVMDT0ervEcujzDxSM3IBqOtqwzeDm9C/aVUxotEXkqSM2oO3Q142mBpgq8Hoc8q7jpr/PJhkZCTx3WvNeG3lBWZsvcKO19vQyMNe44/PU6mZse0qxkYSCwc00vj9hecTiV0Lqthb8r+XGzG/b0NC7j3m/K0Uric+QQYkwNLUmFY1HWnv5UwVe0tdhysUxLsPBKxVpjPq9tB1NMUXdRpMLMHD9z+fsjQzZvVoXwYsO8+4DZfZPqk1tSrbaPTxS47dJDj2EcuHN8e9kvgaL2sisWuRibERzT0daC5Oixqemh3AzAYifzfcxO7pBybmBX7axdaCdWNbMnTlRYasvMiWia2p46KZ5H7xzgOWnbzFYF8PXmxURSP3FIpHFAEThIKYmCvdhiIPglr/Gq08V3qSUvOm5vPrw9R1tWXrpNaoZZkhKy9y837JdnQ9LSolgze3XaGGkzXz+jYs9f2EkhGJXRCepV4fSE+EhCu6jqR47p5RPhaS2EFJ7tsmtUaSYMjKi1y4XfItnreT0xmy8gJ5KpkVI1qIPes6JBK7IDyLVw+QjJRRuyGJOg3mdlClSZFeXsfFlu2TWmNvZcrw1RdZduIWanXxOkndSnrCkJUXyVfJbJ3YWmzX1TGR2AXhWawclRrthpjYq7cD46KPmGtVtmHftPa81Nidrw5HMmb9ZW4lpRd6nUots+1SDK+uuIAsw7ZJIqnrA5HYBeF56vaC+2GGUxTsUSyk3oGaHYt9qY25CUuGNOWzAT5cjkql+7enmLI5kLB7j//T8zdPpebC7Qe8vOwcc38JpXZlG3a83hovV5HU9YGYBBOE5/HuDUc/ghuHodVEXUdTuL/n14uf2EHpGDbcrzq9Grqx9lwUG89H83toIrbmJtRyscHT0YqYBxlEJD4hN1+Nm50FS4Y0pV8Td3GwTo+IxC4Iz+PsBY61lekYQ0jsUafByglcGpTqNk425rzTsx6TOtZmf0g8kYlPuJWUzpWYh1RzsGJM2xo0dLfjhfquYpFUD4l/EUEoTN1ecHkV5DwBcz2eapBlJbHX6ABGmplltbc0ZbhfdY3cSyg7Yo5dEArj3QtUuUpTaH2WegfS7pV4GkYoP0RiF4TCeLb5v/buPUbK6ozj+PdhF0FBFGFbQG5yW1FEihtEqRUVFRDWVJImaI2xbUxabTDxEq1N/2lNmtD0kmLboLFpotSa4KUalUuLtX9IZZWLXBdE1EWRraLIRZDdp3+ctWxg3Z3Zed85M+/8PslkmOXd9zxvIL+cOe95z4FeZ8C2l2JX0rm3/xXec5i/LtmmYBfpSlVPGDsDti8r7adQ334FTh8CA0bHrkQiU7CL5GLcTDjYDLtfj11Jx1pbQ7CPujysTikVTcEukosxM8CqoLFEH1baswEOfQSjrohdiZQABbtILk47C4ZPLd1x9i9v7I7S+Loo2EVyVzsL9m6CfbtiV3KynavC3PXTB8WuREqAgl0kV7Wzw3up9dq/OAzvvKphGPk/BbtIrgaMhoG1YfONUvLuamg5AqOmx65ESoSCXSQftbPCfqKHP4ldyXE7V0GPnjDi0tiVSIlQsIvko3Y2tB6DHStjV3LcW6tg2BToley+pVK+FOwi+RhaB31qSmeN9oMfhamOGl+XdhTsIvnoUQXjroXtK6Dli9jVwNsvh/fRCnY5TsEukq/a2XDk0zDWHttbq8I6NoMnxa5ESoiCXSRfo6ZDde/4wzHuIdjPuSyvbfAk+woKdjP7uZltMLN1ZrbczIYkVZhIyTqlTxjT3vpCCNdY9m6G/U1h022RdgrtsS9094nuPgl4HvhZAjWJlL7xc+HTd+H9tfFqaGx7UErBLicoKNjdfX+7j32AiN0XkSKqnQU9qmHL3+PV0LgcBl8I/QbHq0FKUsFj7Gb2oJm9B9xEJz12M7vNzBrMrKG5ubnQZkXiOu2ssAXd5mfjDMcc+hiaXoOx1xa/bSl5XQa7ma00s40dvK4HcPcH3H0Y8Dhwx1edx90Xu3udu9fV1NQkdwUisZx3fdiO7sNNxW97xz/AW8PUS5ETdBns7j7D3Sd08Hr2hEOXAPPSKVOkBJ07B6xH6LUXW+NLcNpAGDK5+G1LySt0VszYdh/rga2FlSNSRvrWwIhpxR9nb2lb0mDsNdBDM5blZIX+r/hl27DMBuAaYEECNYmUj/H10LwVmrcVr82mNfD5JzBOs2GkY4XOipnXNiwz0d3nuvvupAoTKQvj54b3zUXstW9fFmbkjL6yeG1KWdH3OJFC9BsMwy6Gzc8Ur83GZTD8Euh9RvHalLKiYBcp1Pk3wIcb4cPN6bfVvC08cXrudem3JWVLwS5SqAnzwtDI+r+m39bGpwCD87+dfltSthTsIoXqWwNjroYNT0JrS3rtuMPGpTDym9q0WjqlYBdJwqT5cGAP7Hw5vTb2vAkfbQ/fEEQ6oWAXScK4mdD7TFj/RHptbFwahnzG16fXhmSCgl0kCdW9Qk96y3Pw+f6uj8+XexhfH3UF9BmQ/PklUxTsIkm5cD4cO5zOk6hNDWGZ4Ak3JH9uyRwFu0hShtbBWaPTGY7ZuBSqTtE0R8mJgl0kKWYw6UbY9e9klxho+QI2PR3WhtFDSZIDBbtIki66NeyH+uqi5M655bkw42bSTcmdUzJNwS6SpD4DQq99/d/gwN5kzrn6j9D/nDDzRiQHCnaRpE29HVqOwmsPF36upoawU9LUH2qJXsmZ/qeIJG3gGKidDWsegaOHCjvX6j9ArzM0DCN5UbCLpOHSO+Dwx7B+SffP8elu2PQMTL4ZevVNrjbJPAW7SBqGXxK2rXv1obDjUXeseRhwmHJboqVJ9inYRdJgBpfdFTa7fm1x/r9/eB80/Dnsq9p/RPL1SaYp2EXScu51Ye75P38Bn7yX3+8u/ykc+Qy+dU86tUmmKdhF0mIGs38FOLxwT1jvJRc7X4a1j8GlP4bBE9OsUDJKwS6Spv4jYPr90PhibmvIHD0Ezy0ISxNMvy/9+iSTFOwiaZv6Ixh0Abxwb5jp0plVD8K+XVD/e+h5alHKk+xRsIukraoa6hfB0YPwyFXwwfqTj2lthVcWhlk0F90KI6cVv07JDAW7SDEMmQTfXwZWBY/OgsZlx8fcDzTDYzeEm6wT5sG1D8atVcpedewCRCrG18+HH6yEJd8Jrx7VcGp/OHYkvOb+DibfEm66ihRAwS5STP0Gw60vwrol8NkHYb56y9GwFsygC2JXJxmRSLCb2d3AQqDG3f+bxDlFMqtXX7hYT5NKegoeYzezYcDVwLuFlyMiIoVK4ubpb4B7gRyfvhARkTQVFOxmVg/sdvcO5m+ddOxtZtZgZg3Nzc2FNCsiIp3ocozdzFYCgzr4qweAnwDX5NKQuy8GFgPU1dWpdy8ikpIug93dZ3T0czO7ADgHWG9hetZQ4A0zm+LuexKtUkREctbtWTHu/ibwtS8/m9kuoE6zYkRE4tKTpyIiGZPYA0ruPjKpc4mISPeZ57pGdJKNmjUD7xS94cINBCppqKnSrhd0zZWiXK95hLvXdHVQlGAvV2bW4O51sesolkq7XtA1V4qsX7PG2EVEMkbBLiKSMQr2/HRju/myVmnXC7rmSpHpa9YYu4hIxqjHLiKSMQr2bjCzu83MzWxg7FrSZmYLzWyrmW0ws6fN7MzYNaXFzGaa2TYz22Fm98WuJ21mNszMVpnZFjPbZGYLYtdUDGZWZWZrzez52LWkRcGepwpcf34FMMHdJwKNwP2R60mFmVUBDwGzgPOA+WZ2XtyqUncMuMvdxwNTgdsr4JoBFgBbYheRJgV7/ipq/Xl3X+7ux9o+riYs9pZFU4Ad7r7T3Y8CTwDXR64pVe7+gbu/0fbnzwhhd3bcqtJlZkOB64BHYteSJgV7HvJZfz6jvge8GLuIlJwNvNfucxMZD7n2zGwk8A3gP3ErSd1vCR2z1tiFpEmbWZ8gqfXny0ln1+zuz7Yd8wDhq/vjxaytiKyDn1XEtzIz6wssBe509/2x60mLmc0B9rr762Y2PXY9aVKwn6AS15//qmv+kpndAswBrvLszo9tAoa1+zwUeD9SLUVjZj0Jof64uz8Vu56UTQPqzWw20BvoZ2aPuft3I9eVOM1j76ZKWX/ezGYCvwYud/fM7mloZtWEm8NXAbuBNcCN7r4pamEpstBD+QvwsbvfGbueYmrrsd/t7nNi15IGjbFLVxYBpwMrzGydmf0pdkFpaLtBfAewjHAT8cksh3qbacDNwJVt/7br2nqzUubUYxcRyRj12EVEMkbBLiKSMQp2EZGMUbCLiGSMgl1EJGMU7CIiGaNgFxHJGAW7iEjG/A9zBmmPQOn8xgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x105c74e50>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure();\n",
    "plt.plot(x_test, samples[:,0]);\n",
    "plt.plot(x_test, samples[:,1]);\n",
    "plt.plot(x_test, samples[:,2]);\n",
    "#plt.fill_between(x_test, y_pred - 2*sigma, y_pred + 2*sigma, color=\"#dddddd\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Checking if we can estimate the length scale parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Parameter before learning [2.30258509]\n",
      "Parameters after learning [-0.00784631]\n",
      "Estimated l  [0.99218439] ground truth l  1\n"
     ]
    }
   ],
   "source": [
    "x2_train = x_test\n",
    "y2_train = samples[:,0]\n",
    "\n",
    "kernel2 = RBF(10, (1e-2, 1e2))\n",
    "print \"Parameter before learning\", kernel2.theta\n",
    "gp2 = GaussianProcessRegressor(kernel=kernel2, n_restarts_optimizer=9, alpha=1e-06)\n",
    "gp2 = gp2.fit(x2_train, y2_train)\n",
    "opt_params = gp2.get_params()\n",
    "print \"Parameters after learning\", gp2.kernel_.theta\n",
    "print \"Estimated l \", np.exp(gp2.kernel_.theta), \"ground truth l \", 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Example with multidimensional input"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Parameter before learning [10. 10. 10. 10. 10.]\n",
      "Parameters after learning [1.00000000e+05 3.03896109e+02 2.43607680e+01 1.00000000e+05\n",
      " 4.79395356e-01]\n"
     ]
    }
   ],
   "source": [
    "#Dimensions meaning\n",
    "# 0 - (0/1)\n",
    "# 1 - (0/1)\n",
    "# 2 - int (positive: 1,2,3,4)\n",
    "# 3 - {0,1,2}\n",
    "# 4 - kappa (0.0 - 1.0)\n",
    "Xtrain = np.array([\n",
    "    [0,0,2,0,0.04],  #0\n",
    "    [1,0,2,0,0.04],  #1\n",
    "    [1,1,2,0,0.04],  #2\n",
    "    [1,0,1,0,0.04],  #3\n",
    "    [0,0,1,0,0.04],  #4\n",
    "    [0,0,3,0,0.04],  #5\n",
    "    [1,0,2,2,0.04],  #6\n",
    "    [1,0,2,0,0.00],  #7\n",
    "    [1,0,2,0,0.02],  #8\n",
    "    [1,0,2,0,0.07],  #9\n",
    "    [1,0,2,0,0.04],  # 10\n",
    "    [1,1,2,2,0.04],  # 11\n",
    "    [1,1,2,0,0.0],  # 12\n",
    "    [0,1,1,1,0.1]    # 13\n",
    "])\n",
    "Ytrain = np.array([0.4797,0.4791,0.4868,0.4575,0.4562, #0-4\n",
    "                   0.4649,0.4789,0.4076,0.4730,0.4821, #5-9\n",
    "                  0.4795, 0.4865, 0.4151, 0.4734]) #10-13\n",
    "\n",
    "\n",
    "l = np.full(5,10.0)\n",
    "# se_kernel = RBF(l, (1e-5,1e2))\n",
    "se_kernel = Matern(l)\n",
    "#se_kernel = Matern(l, nu=1.2)\n",
    "print \"Parameter before learning\", np.exp(se_kernel.theta)\n",
    "\n",
    "\n",
    "## learning part\n",
    "gpReg = GaussianProcessRegressor(kernel=se_kernel, n_restarts_optimizer=9, alpha=1e-06)\n",
    "gpReg = gpReg.fit(Xtrain, Ytrain)\n",
    "print \"Parameters after learning\", np.exp(gpReg.kernel_.theta)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.37112173 0.37299963]\n"
     ]
    }
   ],
   "source": [
    "Xtest = np.array([[1,1,2,0,0.01], #10\n",
    "                  [1,1,2,2, 0.04] #11\n",
    "                 ])\n",
    "\n",
    "# expected output [0.358, 0.356]\n",
    "y_pred, sigma = gpReg.predict(Xtest, return_std=True)\n",
    "print y_pred"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Model estimation\n",
    "Finding out which kernels to use based on Leave One Out Cross validation.\n",
    "For that we construct the set of potential kernels to choose from and set of training sets based on the number of training data point we have.\n",
    "Then for every kernel and for every set we estimate the -log marginal likelihood, which is also a generalization error"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "we have  14  training points\n"
     ]
    }
   ],
   "source": [
    "# creating kernels set\n",
    "l = np.full(5,10.0)\n",
    "kernels = []\n",
    "kernels.append(RBF(l, (1e-5,1e2)))\n",
    "kernels.append(Matern(l))\n",
    "kernels.append(Matern(l, nu=0.5))\n",
    "\n",
    "from itertools import combinations\n",
    "N = Xtrain.shape[0]\n",
    "print \"we have \", N, \" training points\"\n",
    "indices = range(N)\n",
    "sets = list(combinations(indices, N-1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[array([344.17857742]), array([-29.85014166]), array([-53.71425435])]\n",
      "Kernel with min error is  2  with error  -53.71425434688109\n"
     ]
    }
   ],
   "source": [
    "from scipy.stats import norm\n",
    "\n",
    "gen_error = [] # empirical generalization errors or empirical risks\n",
    "for kernel in kernels:\n",
    "    gen_error_per_kernel = 0 # empirical generalization error for a kernel\n",
    "    for i,s in enumerate(sets):\n",
    "        # selecting those training points that are in the set s\n",
    "        X = Xtrain[s,:]\n",
    "        Y = Ytrain[[s]]\n",
    "        ## finding missing point\n",
    "        missing_idx = sum(range(N)) - sum(s)\n",
    "        # set Xl - X left out\n",
    "        Xl = [Xtrain[missing_idx]]\n",
    "        Yl = Ytrain[missing_idx]\n",
    "        \n",
    "        ## fit the gaussian process with a given kernel\n",
    "        gp = GaussianProcessRegressor(kernel=kernel, n_restarts_optimizer=9, alpha=1e-06)\n",
    "        gp = gp.fit(X, Y)\n",
    "        y_pred, sigma = gp.predict(Xl, return_std=True)\n",
    "#         loss = -np.log(norm.pdf(Yl, y_pred, sigma))\n",
    "        loss = -norm.logpdf(Yl, y_pred, sigma)\n",
    "        # print loss\n",
    "        gen_error_per_kernel = gen_error_per_kernel + loss\n",
    "    gen_error.append(gen_error_per_kernel)\n",
    "\n",
    "# generrpr /= N_test\n",
    "# finding the kernel with minimal average generalization error\n",
    "print gen_error\n",
    "min_error =  np.min(gen_error)\n",
    "min_error_kernel = np.argmin(gen_error)\n",
    "print \"Kernel with min error is \", min_error_kernel, \" with error \", min_error "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Searching for input points that maximize the value of the output\n",
    "For now we know what kernel fits our data the best + we know the parameters for that kernle to use.\n",
    "Now we can estimate the posterior mean function $\\mu_\\mathrm{post}$ and covariance for the GP. The posterior mean function gives the outputs of the model for every possible input point. So, if now we need to estimate what the value (y-coordinate) for any other input point x is, we should just compute $y = \\mu_\\mathrm{post}(x)$. And it will work, because we already estimated the model given our training points.\n",
    "\n",
    "\n",
    "Now, let's assume we want to know what should the input point $x^*$ be so that our $\\mu_\\mathrm{post}$ has the biggest value (because we need the biggest value for some other purposes). For that we need to search for the maximum of the $\\mu_\\mathrm{post}$. This is a naive approach to do that. But what the GP posterior also tells us is not just the function $\\mu_\\mathrm{post}$, but also the confidence about this function. \n",
    "This is important to have, because the mean function estimate can have a maximum peak somewhere, but we are very uncertaint about that peak. This is not a desirable situation. We would prefer a slightly lower peak, but with more confidence. \\todo{put a picture here}\n",
    "Luckily, there exists a concept called **acquisition function**. This function combines the mean estimates and the uncertainty about that estimates. There are a lot of different ways to combined that information. The most popular one is called **expected improvement** function and luckily if we use Gaussian Processes it has a close form solution.\n",
    "see here for details: https://thuijskens.github.io/2016/12/29/bayesian-optimisation/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(1, 5)\n",
      "Expected improvement [-0.00818707]\n",
      "posterior mean  [0.47710448]  sigma [0.03117798]\n",
      "params  [11.51292546 10.4491895   7.98426886 11.51292546  3.90194279]\n",
      " Gaussian 0.7090134490872383 \n",
      " camAware 0.6893486565465319 \n",
      " kernelSize  3.0 \n",
      " vig_opt  0.6637448026125269 \n",
      " kappa 0.1\n"
     ]
    }
   ],
   "source": [
    "## Functions taken from https://github.com/thuijskens/bayesian-optimization/blob/master/python/gp.py\n",
    "from scipy.stats import norm\n",
    "from scipy.optimize import minimize\n",
    "\n",
    "def expected_improvement(x, gaussian_process, evaluated_loss, greater_is_better=False, n_params=1):\n",
    "    \"\"\" expected_improvement\n",
    "    Expected improvement acquisition function.\n",
    "    Arguments:\n",
    "    ----------\n",
    "        x: array-like, shape = [n_samples, n_hyperparams]\n",
    "            The point for which the expected improvement needs to be computed.\n",
    "        gaussian_process: GaussianProcessRegressor object.\n",
    "            Gaussian process trained on previously evaluated hyperparameters.\n",
    "        evaluated_loss: Numpy array.\n",
    "            Numpy array that contains the values off the loss function for the previously\n",
    "            evaluated hyperparameters.\n",
    "        greater_is_better: Boolean.\n",
    "            Boolean flag that indicates whether the loss function is to be maximised or minimised.\n",
    "        n_params: int.\n",
    "            Dimension of the hyperparameter space.\n",
    "    \"\"\"\n",
    "\n",
    "    x_to_predict = x.reshape(-1, n_params)\n",
    "\n",
    "    mu, sigma = gaussian_process.predict(x_to_predict, return_std=True)\n",
    "\n",
    "    if greater_is_better:\n",
    "        loss_optimum = np.max(evaluated_loss)\n",
    "    else:\n",
    "        loss_optimum = np.min(evaluated_loss)\n",
    "\n",
    "    scaling_factor = (-1) ** (not greater_is_better)\n",
    "\n",
    "    # In case sigma equals zero\n",
    "    with np.errstate(divide='ignore'):\n",
    "        Z = scaling_factor * (mu - loss_optimum) / sigma\n",
    "        expected_improvement = scaling_factor * (mu - loss_optimum) * norm.cdf(Z) + sigma * norm.pdf(Z)\n",
    "        expected_improvement[sigma == 0.0] == 0.0\n",
    "\n",
    "    return -1 * expected_improvement\n",
    "\n",
    "def get_next_param_set(acquisition_func, gaussian_process, evaluated_loss, bounds, greater_is_better=False, n_restarts=25):\n",
    "    \"\"\" sample_next_hyperparameter\n",
    "    Proposes the next hyperparameter to sample the loss function for.\n",
    "    Arguments:\n",
    "    ----------\n",
    "        acquisition_func: function.\n",
    "            Acquisition function to optimise.\n",
    "        gaussian_process: GaussianProcessRegressor object.\n",
    "            Gaussian process trained on previously evaluated hyperparameters.\n",
    "        evaluated_loss: array-like, shape = [n_obs,]\n",
    "            Numpy array that contains the values off the loss function for the previously\n",
    "            evaluated hyperparameters.\n",
    "        greater_is_better: Boolean.\n",
    "            Boolean flag that indicates whether the loss function is to be maximised or minimised.\n",
    "        bounds: Tuple.\n",
    "            Bounds for the L-BFGS optimiser.\n",
    "        n_restarts: integer.\n",
    "            Number of times to run the minimiser with different starting points.\n",
    "    \"\"\"\n",
    "    best_x = None\n",
    "    best_acquisition_value = 1\n",
    "    n_params = bounds.shape[0]\n",
    "\n",
    "    for starting_point in np.random.uniform(bounds[:, 0], bounds[:, 1], size=(n_restarts, n_params)):\n",
    "\n",
    "        res = minimize(fun=acquisition_func,\n",
    "                       x0=starting_point.reshape(1, -1),\n",
    "                       bounds=bounds,\n",
    "                       method='L-BFGS-B',\n",
    "                       args=(gaussian_process, evaluated_loss, greater_is_better, n_params))\n",
    "\n",
    "        if res.fun < best_acquisition_value:\n",
    "            best_acquisition_value = res.fun\n",
    "            best_x = res.x\n",
    "\n",
    "    return best_x\n",
    "\n",
    "\n",
    "# iteratively search for the improvement\n",
    "# fit the model, which previously performed the best\n",
    "kernel = kernels[min_error_kernel]\n",
    "model = GaussianProcessRegressor(kernel=kernel, n_restarts_optimizer=25, alpha=1e-06)\n",
    "model.fit(Xtrain, Ytrain)\n",
    "bounds = np.array([[0, 1],\n",
    "                  [0, 1],\n",
    "                  [1, 3],\n",
    "                  [0, 2],\n",
    "                  [0.0, 0.1]]\n",
    "                 ) # bounds for the optimizer\n",
    "next_sample = get_next_param_set(expected_improvement, model, Ytrain, bounds, greater_is_better=True, n_restarts=100)\n",
    "next_sample = next_sample.reshape(-1,1)\n",
    "next_sample = np.transpose(next_sample)\n",
    "print next_sample.shape\n",
    "\n",
    "impr = expected_improvement(next_sample, model, Ytrain, greater_is_better=True, n_params=5)\n",
    "print \"Expected improvement\", impr\n",
    "y_pred, sigma = model.predict(next_sample, return_std=True)\n",
    "print \"posterior mean \", y_pred, \" sigma\", sigma\n",
    "print \"params \", model.kernel_.theta\n",
    "\n",
    "print \" Gaussian\",next_sample[0,0], \"\\n camAware\",next_sample[0,1], \"\\n kernelSize \",next_sample[0,2], \"\\n vig_opt \", next_sample[0,3], \"\\n kappa\", next_sample[0,4] "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
